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We study two-phase stratified flow where the bottom layer is a thin laminar liquid and 
the upper layer is a fully-developed gas flow. The gas flow can be laminar or turbu- 
lent. To determine the boundary between convective and absolute instability, we use 
Orr-Sommcrfcld stability theory, and a combination of linear modal analysis and ray 
analysis. For turbulent gas flow, and for the density ratio r = 1000, we find large regions 
of parameter space that produce absolute instability. These parameter regimes involve 
viscosity ratios of direct relevance to oil/gas flows. If, instead, the gas layer is laminar, ab- 
solute instability persists for the density ratio r = 1000, although the convective/absolute 
stability boundary occurs at a viscosity ratio that is an order of magnitude smaller than 
in the turbulent case. Two further unstable temporal modes exist in both the laminar 
and the turbulent cases, one of which can exclude absolute instability. We compare our 
results with an experimentally-determined flow-regime map, and discuss the potential 
application of the present method to non-linear analyses. 



1. Introduction 



We investigate linear absolute and convective instability for a liquid film sheared by 
laminar and turbulent gas streams in a channel. In the oil/gas industries, this approach 



serves as a model that can be used to predict the onset of droplet entrainment (Hall- 
Taylor & Hewitt 19701). The motivation for this work is twofold: previous work on the 



Boomkamp & Miesen 



1996; 



Boomkamp et at 



1997 



r analysis ( 


Miesen & Boersma 1995 


Ozgen et al. 


1998 


), while previous 



that are relevant to the oil-and-gas industries, and which are found herein to be absolutely 
unstable. The current work aims to fill in these two gaps in the literature and introduces 
modifications and extensions of existing methodologies (developed previously for single- 
phase flows or for temporal stability analysis only) that are potentially of interest in 
other areas. 

The route to droplet entrainment from a liquid layer into a gas stream in pipe and 
channel flows is still unclear. The idealized system of fully-developed flow with a flat 
gas/liquid interface is linearly unstable to infinitesimally small perturbations. For a lam- 
inar base state, a stratification in dynamic viscosity produces instability; for a turbulent 



base state, the mechanism of Miles (1957) may also dominate for deep liquid layers at 
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large Froude numbers (O Naraigh et al. 2011&). Other mechanisms, such as a Tollmicn 



Schlichting mode in the liquid (Boomkamp & Miesen 1996) may also be important in 



certain parameter regimes. Although linear (temporal) instability is arguably a necessary 
condition for droplet entrainment in the gas, it does not provide much insight regard- 
ing whether a localized disturbance grows whilst it is merely convected downstream or 
whether instead the disturbance destabilizes the entire system. To answer these ques- 
tions a spatio-temporal analysis is required. A recent spatio-temporal analysis by | Vallu"rT| 
et al. (2010) has revealed a region in parameter space wherein the laminar base state is 



absolutely unstable, indicating that the system does not merely act as an amplifier but 
also as a generator of disturbances. This was found to include a large range of practically 
useful viscosity ratios, but was limited to a density ratio of O(l). 

The laminar density-matched problem studied by Valluri et al.\ (2010) has only lim- 
ited applicability in oil/gas transport, where the density ratio is large, and where the 
operating conditions produce turbulence in the gas layer, or in both layers (see, e.g., 



the visualisation of droplet entrainment events by Lecoeur et al. (2010)). The present 



study therefore investigates the corresponding problem for a turbulent base state, but 
the l aminar case is also revisited. Although replaceme nt of the base state by a turbulent 



one ( |6 Naraigh & Spelt|2010( |0 Naraigh et al. |2011 al in a linear modal spatio-temporal 



analysis may seem trivial, the results turn out to be difficult to interpret, due to the pres- 
ence of multiple unstable modes. Therefore, in this study, we have developed a twin-track 
approach, in which modal analysis and ray analysis are combined to locate and charac- 



terize absolute instability. The ray analysis used herein extends the work of Delbende 



et al. ( 1998[ ) for single-phase flows. This approach yields surprising results. In particular, 
it has revealed significant regions in parameter space where the turbulent base state is 
absolutely unstable for large density ratios. This also holds for the laminar base state, 



thereby contradicting Valluri et al. (20101, who only found absolute instability for den- 



sity ratios of 0(1). We have therefore revisited the spatio-temporal work of Valluri et al. 



(2010), and have established using both ray and modal analyses that the laminar system 



is indeed absolutely unstable at large density ratios for a substantial range of viscosity 
ratios and liquid-film depths. Reasons for the oversight in the previous work are given. 

The paper is organized as follows. The turbulence model and the linear stability anal- 
ysis are formulated in Section [2] We discuss some theoretical and numerical aspects of 
the linear stability analysis in Section [3] paying close attention to the development of a 
ray analysis for two-phase flows. We apply this technique to the turbulent base state in 
Section [4j while the laminar case is revisited in Section [5] In Section [6] we argue for the 
importance of using the ray analysis and the modal analysis simultaneously, for complete 
and accurate results. We also compare the flow-regime boundaries identified herein with 
those found in experiments, and discuss the generalisation of our work to non-linear and 
non-parallel flows. 



2. Linear stability analysis 

In thi s section we review a model of turbulent channe l flow used elsewhere by the 



authors (0 Naraigh & Spelt 2010 Naraigh et q/.||2011«). This is a Reynolds-averaged 



model describing co-current flow in a stably-stratified system where the upper layer 
is a turbulent gas and the lower layer a laminar liquid film. We also recall the Orr- 
Sommerfeld technique used to determine the stability of the interface in this two-layer 
system. For reference, typical values of the problem parameters are given in Table [I] 
where the subscripts G, L indicate the gas and liquid, respectively. 
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Physical Quantity 


Symbol 


Numerical value 


Gas-layer dynamic viscosity 


Hg 


1.8 x 10" 5 Pa- s 


Liquid-layer dynamic viscosity 


PL 




Viscosity ratio 


m = pl/pg 


50-10000 


Gas-layer density 


PG 


lkgm~ 3 


Liquid-layer density 


PL 




Density ratio 


r = phj pa 


10 3 


Liquid film thickness 


d L 


10" 3 -10 -2 m 


Gas-layer depth 


d G 




Ratio of layer depths 


e = d L /d G 


0.01-0.2 


Acceleration due to gravity 


9 


9.8ms" 2 


Surface tension 


7 


0.07 Nm" 1 



Table 1. Table of parameters in the two-phase problem and their typical values. The range 
of values of the liquid-layer dynamic viscosity pl, the liquid-layer density pL, and the liquid- 
film thickness dL can be backed out from the gas-layer analogues and the ratios m, r, and e, 
respectively. 



2.1. The base state 

We consider a flat-interface base state in two-layer stratified flow (Figure [T]). The bottom 
layer is a thin, laminar, liquid layer, and the top layer is gaseous, turbulent and fully- 
developed. A pressure gradient is applied along the channel. The base-state profile of 
the system is a uni-directional flow in the horizontal, x- direction as a function of the 
cross-flow coordinate z. In the bottom layer, the profile is determined by balancing the 
viscous and the pressure forces; in the top layer, the viscosity in the balance law must 
be supplemented by the turbulent eddy viscosity: 

dU 



dz 



to = n 



dP 

dL 



z, 



(2.1) 



where Uo(z) is the base-state velocity in the gas. 
dP/dL is the applied pressure gradient. Moreover, 
stress due to the averaged effect of the turbulent fluctuating velocities, u and w . In 



T\ is the interfacial shear stress, and 
To = —pg(u'w') is the turbulent shear 



channel flows, it is appropriate to model this term using an eddy- viscosity model ( Monin 



fc Yaglom 1971). In mixing- length theory, the eddy viscosity depends on the local rate 
of strain ( Bradshaw||1974 ), which means that the t urbulent shear stress dep ends on the 



square of the rate of strain. Instead, as in the work of O Naraigh et al. ( 2011a 



and 



Biberg 



(2007), we use an interpolation function for the eddy viscosity. This mimics the ordinary 
mixing-length theory near the interface and near the wall, and transitions smoothly from 
the near-wall and near-interfacial regions to the zone surrounding the gas centreline. 
Thus, the turbulent shear stress is linear in the rate of strain, and 

dU 
dz ' 

where fix is the eddy viscosity and k is the von Karman constant, taken as 0.4. Addition- 
ally, U* w is the friction velocity at the upper wall. The corresponding stress is t w , with 
U* w = \j \t w \/pg- Similarly, the interfacial friction velocity is defined as = Pg - 
The function T is the interpolation function described in O Naraigh & Spelt (20101 



T = Mt- 



P-t = npcdaU^J 7 (s) rfc (s) ip w (1 - s) . 



z/d c 



(2.2) 
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Liquid 



-z=0 



Figure 1. A schematic diagram of the base flow. The liquid layer is laminar, while the gas layer 
exhibits fully-developed turbulence, described here by a Reynolds-averaged velocity profile. A 
pressure gradient in the a;-direction drives the flow. 



and 



Q Naraigh et al. (2011a 



F(s) = s(l-s) 



\R\ 5 / 2 (1 - sf 



R 2 (l-sf + Rs(l-s) + s 2 



R = t;/t w 



(2.3) 



Finally, tpi and ^ w are interface and wall functions respectively, which damp the effects 
of turbulence to zero rapidly near the interface and the wall. These are given below. 

We non-dimcnsionalize the problem on the gas-layer depth da, the gas-layer density 
pG, the gas-layer viscosity po, and the velocity scale U p , where 

PG U 2 = d G \dP/dL\. 



Then, integration of Equation (2.1 1 yields the non-dimensional base state: 

(z + e) 



i'L 



U (z) 



1 1 2 n , R £ l \ , R £ l r z 



ffl(«)V>i(*)V'w(l-*) : 



(2.4) 



where e = d^/dc, Re = pcUpdc/ Pg, an d where Re* = (U*/U p )Re. Knowledge of Re* 
amounts to knowledge of the interfacial shear stress. This is not known a priori as a 
function of the externally-imposed parameters. However, it is available within the model, 
and the root-finding procedure 



U (1; Re*) = 



(2.5) 



yields Re* as a function of the parameters (Re, e, pl/pg)- For completeness, we also list 
the interfacial and wall functions: 



A (s) 



(2.6a) 
(2.66) 



where Ca is a constant fixed such that the interfacial and wall viscous sublayers are five 
wall-units in extent (Pope 2000). The functional forms for T and the wall functions are 
confirmed by the excellent agreement bet ween the model predictio ns of the base state 
and experiments and numerical simulation ( O Naraigh et al. |2011 al). Having constituted 
the base state, we now introduce the theory necessary to determine its stability. 
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2.2. The perturbation equations 

We base the dynamical equations for the interfacial motion on the Reynolds-averaged 
Navier-Stokes (RANS) equations. The turbulent velocity is decomposed into averaged 
and fluctuating parts. We make the quasi-laminar approximation, which means that 
the fluctuations are only considered in the base state, where they are modelled using 
the eddy viscosity (Section 2.1). Before deriving equations for perturbations induced by 
small waves at the interface, we discuss the dynamics of these perturbations with respect 
to their interactions with the turbulent eddies in the flow. 

In a realisation of the three-dimensional turbulent two-layer flow with small-amplitude 
waves, a Fourier mode decomposition can be made of the interface elevation and field 
variables. Here, ensemble-averaged Fourier modes are assumed to be predominantly two- 
dimensional. The correspondingly averaged equations of motion are linearized in terms 
of wave amplitude. The linearized problem contains wave-induced Reynolds stress terms 
(WIRSs) , but these have b een found recently not to be significant in two-layer flows 
such as those studied here (O Naraigh et al. 2011a 6). This can also be seen from an 



order-of-magnitude estimate of the WIRSs terms compared to inertial terms in the per- 
turbation momentum equation. Further theoretical justification exists for the case of 
viscosity-contrast instability, where the instability is dominated by conditions close to 
the interface, a zone where the perturbation turbulent stresses are damped rapidly to 
zero by the existence of viscous sublayers. In practical terms, the quasi-laminar approx- 
imation, wherein the WIRSs are ignored and the effect of turbulence is assumed to be 
entirely through the base-state velocity profile, while brutal in its simplicity, yields simi- 
lar results to other turbulence models that explicitly include the WIRSs. It also predicts 
critical Reynolds numbers for the onset of wavy flow that agree with the laboratory ex- 
periments^ Cohfin_^FIairraJ|Y J 1965 ) and Craik (1966). The reader is referred to the 

|2011a|&| for further details. 



papers by O Naraigh et al. 



Here, we study linear spatio-temporal instability. Although the above-mentioned prior 
findings are limited to temporal Fourier modes, there is an analy tical connection between 
spatio-temporal and temporal modes (6 Naraigh et al. 2012) (see also Appendix |a| , 
such that the properties of the temporal study are inherited by the spatio-temporal one. 
Therefore, we make here also the quasi-laminar approximation. Thus, a small disturbance 
z = rj{x, t) centred around the flat interface z — gives rise to disturbances in the velocity 
and pressure fields, which satisfy the following linearized equations of motion in the j th 
phase (j — L,G): 



d r tt 9 r 

g- t Su + U -Su- 



d 



dx 



dz 

rr 9 r ' 

- U — dw 
ox 

5u + -^—Sw = 0, 
oz 



d_ 

dx 
d_ 
dz 



Sp - 
Sp - 



Re 
mj_ 
Re 



cf_ 
dx 2 

dx 2 



dz 2 
dz 2 



Su, 
Sw, 



(2.7a) 
(2-76) 
(2.7c) 



where (rL,rc) = {r, 1) and (mi,, me) — (m, 1). Using the incompressibility condition, 
this system of equations reduces to a single equation in the streamfunction <j). Further 
simplification occurs when the streamfunction is written as a sum of normal modes: 



4>(x,z,t) 



1 

2^ 



dae iax cf> a (z,t), 



(2.8) 



(2.9) 
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which in turn can be written in Laplace-transform notation: 

where is the Bromwich contour ( Arfken fc Weber]|2001 ). If the w-singularities in the 
function 4> au (z) lie below the real axis in the complex-w plane, then the integral (2.9 1 
is an ordinary double Fourier integral. Using the Fourier and Laplace decompositions, 
Equations (2.7 1 reduce to the Orr-Sommerfeld equation: 

d 2 U 



Re 



(2.10) 



where D = d/dz. Equation ( |2.10[ ) only holds in the interior parts of the domain, z € 
(— e, _ ) U (0 + ,l). To close the Equation (2.10), no-slip and no-penetration conditions 
are applied at z = — e and z = 1: 



bauj (~e) = U(j) auJ (-e) = (pauj (1) = D0 QW (1) = 0, 



(2.11) 



and the streamfunction is matched across the interface z = 0, where the following con- 
ditions hold (we use the notation c = ui/a): 



m (D + a 2 



yL = <PG, 
- (D 2 H 



<t>G 

c-U 

, 2\ 



(dU 




dU 




\ dz 


0+ 


dz 


J 



o?)<t> 



G ■ 



(2.12a) 
(2.126) 
(2.12c) 



77i (D 3 (f) L - 3a 2 D0 L ) + \arRe (c - U ) D0 L + iarRe 



dUo 

dz 



> L -^(F + * 2 S)^ L 



(D 3 <j> G - 3a 2 D0 G ) + iaRe (c - U ) D0 G + iaRe 



c-U 

dUo 

dz 



b G . (2.12d) 



o+ 



Here F and S denote parameters that encode the effects of gravity and surface tension, 
respectively; they are defined here for the first time as 



F = 



S 



7 1 



2 ^ := F (r-l)/Re 2 



I^gIpgAg Re 2 



So/Re 2 



(2.13) 
(2.14) 



where g is acceleration due to gravity and 7 is surface tension. The appropriate range of 
values for F and So is discussed in Sec. [4] We abbreviate the Orr-Sommerfeld (or OS) 
equation (2.10) and the matching conditions ( |2.11 )— ( 2.12 ) using operator notation, 



C 



= \uM aL 



(2.15) 



This equation amounts to an eigenvalue equation, which we solve numerically by intro- 
ducing a trial solution: 



J OL UJ (z) « 2J a n T n 



n=0 
N 2 



2z 
e 



1 



-e^z^O, 



b aul {z)^^b n T n {2z-l), 0<2<1, 



(2.16a) 
(2.166) 
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where T n {-) is the n th Chebyshev polynomial. We substitute Equation (2.16) into Equa- 
tion (2.15) and evaluate the result at N 2 — Q interior points. The ansatz (2.16) is also 
substituted into the eight boundary and interfacial conditions. This yields N\ + N% +2 
linear equations in as many unknowns. In matrix terms, we have to solve 



VjjM„ 



(2.17) 



where v = (ao, .., , bo, .., 6jv 2 ) . Such an equation is readily solved using linear-algebra 
packages. This method is described in more detail and its implementation is tested against 



benchmarks in another paper by the present authors (Valluri et al. 2010). The number 



of collocation points (N± + 1,N% + 1) is adjusted until convergence is achieved. The 
application of this numerical method will be the subject of the following sections. 



3. Further numerical methods and postprocessing 

In this section we revisit the basic definition of absolute instability, namely that the 
streamfunction response to a localized disturbance should grow exponentially in time at 
the origin of the disturbance. Solving the associated Cauchy problem gives a quick and 
clear method to characterize the instability. This approach also enables us to pinpoint 
the source of the instability through an energy-budget analysis. We also review herein 
an equivalent method to determine absolute instability, namely modal analysis. 

3.1. Modal analysis 



A purely temporal analysis involves the solution of the eigenvalue problem (2.15), where 
we write a = a x + iaj, u) = ui r + icjj, for a = a r only. This gives a dispersion relation 

(w[ omp (a r ),w[ cmp (a r )) _ ( Wr ( ar)Q! . = 0),Wi(a r ,ai = 0)) , 

with associated group velocity c g = dw* omp /da r . The pair (a r ,o;; emp (a r )) that maximizes 
u)\ is called the most dangerous mode. The flow is linearly unstable if Wj Cmp > for 
the most dangerous mode. Unstable parallel flows are further classified as convectively 
unstable if initially localized pulses are amplified in at least one moving frame of reference 
but are damped in the laboratory frame, and absolutely unstable if such pulses lead to 
growing disturbances in the entire domain in the laboratory frame. To describe such an 
instability, we use the description of Huerre & Monkewitz ( 1990[ ). An unstable parallel 
flow is absolutely unstable if the following criteria have all been met: (i) wio := Wj(ao) > 0, 
where ao is the wave number at which the complex derivative dio/da is zero, (ii) the 
corresponding saddle point ao m t ne complex a-plane is the result of the coalescence of 
spatial branches that originate from opposite half-planes at a larger and positive value 
of u)i and (iii) the saddle point pinches at w;o; this is verified by locating a cusp at cjjo in 
the complex w plane ( Lingwood|1997 Valluri et aL|2010| and ensuring that the complex 
wave number corresponding to the pinching point coincides with ao- 



3.2. Ray analysis 



The linear stability equations (2.7) in streamfunction form (Equation (2.10)), together 



with the boundary and initial conditions ( 2.11 )-( 2.12) can be neatly encoded in linear- 
operator form: 

£4> = Md t 4>, (3.1) 
where we study <p(x,z,t), the filtered streamfunction containing only positive real wave 



numbers (Huerre & Monkewitz 1990); here C and M. are linear operators (cf. Equa- 
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tion (2.15)). When an impulsive, localized force is applied to the streamfunction, Equa- 
tion (3.1) is modified: 

C4> + S(x)S(z)6(t) = Mdtfa (3.2) 
in terms of the Dirac delta function S(.). The solution of Equation (3.2) (determined 



here in a domain with periodic boundaries at x = ±L x /2) can be used to characterize 
the instability that develops from the impulse. We use the following algorithm developed 
by Delbende & Chomaz (19981 and Delbende et al. (1998), here applied to two-phase 



flows: 



(a) Compute the complex-valued filtered streamfunction <f> (x, z, i); 
(6) Form the L 2 -norm 



A(x,t) 




(3.3) 



(c) Examine the norm along rays, A(v,t) = A(x = vt,t). If A(0,t) is a decreasing 
function of time, the instability is convective. 

It suffices to consider positive and zero ray velocities only, since this enables a classifica- 
tion of the instability as either absolute or convective and, furthermore, in the convective 
case, gives information about the leading- and trailing-edge velocities of the downstream- 
propagating disturbance. 

Additional information can be extracted from the evolution of the norm, provided 
the contributions to the growth of the streamfunction are dominated by a single eigen- 



mode. This caveat does not appear in the single-phase work of Delbende &: Chomaz 



(1998) and Delbende et al. (19981: those problems contain a simple means of projecting 



the streamfunction on to the eigenmode of interest; the spatial symmetries that pro- 
duce this projection do not exist in the current problem, and this approach is therefore 
not applicable. We therefore assume that the evolution is dominated by a single eigen- 
mode, and justify this assumption a posteriori. Thus, along rays x = vt, we assume that 
A(x = vt,t) ~ t-i/V^*, where a is the spatiotemporal growth rate of the dominant 
eigenmode. Therefore, we extract the finite-time estimate of the spatiotemporal growth 
rate as follows: 



a(v) = 



hi A (vt2, t%) — In A {vt\, t±) 1 In £2 — hiii 



t 2 ~ti 



t 2 -h 



(3.4a) 



where t\ and ti are large but finite times and ti > t\. The complex wave number and 
frequency along the ray x = vt also follow from this analysis ( Delbende fc Chomaz|1998~ 
Delbende efaT|l998| ): 



Otr(v) = 5ft 



da 

dv ' 
a(v) + ait?, 

6 dx 



(3.46) 
(3.4c) 

(3.4d) 



where Equation (3.4d) holds because the right-hand side is independent of time as t — > 00. 



3.3. Transient direct numerical simulations 



Transient direct numerical simulation (DNS) of Equation (3.2) is complicated by the fact 



that the operator M. is non-invertible. To solve this equation in an optimal way, we have 
developed our own numerical method, which we outline here. As in Section [2j we write 
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4> (z, t) as a finite sum of Chebyshev polynomials: 



or more compactly, 



(3.5) 



a>0 n 



We substitute Equation (3.5) into Equation (3.2). This yields the following equation for 
the normal mode a: 

di? 

M a -^=L a v a , t>0, (3.6) 



where M a and L a are the Orr-Sommerfeld matrices described in Equation ( 2.17[ ). The 
matrix M a is not invertible: it has rows of zeros corresponding to the no-slip boundary 
conditions, the continuity of the streamfunction at the interface, and the continuity of 



the tangential stress at the interface. Equation (3.6) is therefore a differential algebraic 
equation (DAE). 

There are several standard methods for solving DAEs with computational software 
packages ( Shampine et aL| [ l9 99). For a singular matrix M, the DAE M(t,y)y' = f(t,y) 
has a solution only when the initial condition y Q is consistent, that is, if there is an 
initial slope y p o such that M(t ,y a )y p0 = f(t ,y ). In general, computational packages 
for solving DAEs demand not only that the initial data be consistent, but also that 
the slope be prescribed as an input to the numerical solver ( Shampine et aL||1999 |. We 
develop herein a numerical method for linear DAEs that removes the necessity to specify 
the slope. Moreover, long-time integrations of DAEs using computational packages can 
be costly, especially for the modal decomposition (3.5), which contains a large number 



of wave numbers. Thus, we resort to a semi-analytical solution method that holds for 



constant-coefficient DAEs such as Equation (3.6). We re- write Equation (3.6) as 

d 



dt 



M a v a — L a v ai 



and integrate it over a numerical time step At: 
M a v a (t + At) - M a v a (t) = 



t+At 



L a v a 



ids. 



For a sufficiently small timestep, the integral on the right-hand side may be approximated 
by the trapezoidal rule. Thus, the equation is approximated by 

M a v a (t + At) - M a v a (t) = \At [L a v a (t + At) + L a v a (t)] . 

Re-arrangement gives 

[M a - \AtL a ] < +1 = [M a + \AtL a ] <, 

where v^ +1 — v a (t + At) and = v a (t). The solution at an arbitrary time p is therefore 
available from the initial condition using only three matrix operations: 

-l 



\M n 



\AtL 



AtL a 



(3.7) 



Equation (3.7) becomes an exact solution to the DAE ( |3.6| in the limit when At — > 0: 
v a (t) = X a (t) v 0a , X a (t) = lim XI. 



At-j-0 
t=pAt 
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.(ST) 



.(Modal) a(ST) a(Modal) 



400 
1300 



1.52 
1.39 



1.52 
1.39 



20 
17 



21 

17 



Table 2. Comparison between modal and spatio-temporal analyses at various values of m , with 
r = 1000, e = 0.1, and Re — 2500. The parameters S and F are given in Equation (4.1 1. The 
DNS parameters are Ni = 21, N 2 = 60, L x = 150, Ax = 0.01, and At = 0.01. 



Finally, to approximate a delta-function impulse, we fix the initial condition v 0a as 
follows: 



>(x,z,t = 0) = ^V^e" 



a>0 




where N n = 1/ir if n = and 2/-7T otherwise. The coefficients w x and are the widths 
of the approximate delta functions in the x- and z-directions respectively. 

This solution method is appropriate for the kind of long-time simulations performed 



herein because the time required on a computer to implement Equation (3.7 1 depends 



only very weakly on p: the computation time required to raise an arbitrary squa re matrix 
to the p th power depends (in a worst-case-scenario) only logarithmically on p ( Ar & Cai 



1994). Thus a DNS of 1000 time units takes (at most) only three times longer than a 



DNS of 10 time units. 

We have validated the solution method by computing the growth rate for delocal- 
ized monochromatic initial conditions, and checking that the dispersion curve agrees 
with standard (temporal) eigenvalue analysis: the results are identical in each case and 
are not reported further here. We have further validated the spatio-temporal transient 



DNS by using Equations (3.4) and computing the growth rate and wave number at the 



most dangerous spatio-temporal mode. This is necessarily the most dangerous tempo- 



ral mode (Huerre & Monkewitz 1990). Thus, if our transient DNS is correct, then the 



maximum growth rate computed in this manner must agree with standard temporal 
eigenvalue analysis. The results for this comparison are shown in Table [2j and confirm 
the correctness of our transient DNS. We now apply this method to classify completely 
the instability in turbulent two-phase stratified flow. 



4. Results for the turbulent base state 

4.1. Modal analysis and flow-regime maps 

In this section we examine the spatio-temporal instability wherein the upper layer is 
turbulent. We base the parameter values on an upper layer of air of depth da = 5 cm, 
together with the values of surface te nsion and gravitational a cceleration given in Table [T] 
(these parameter values were used in O Naraigh et al. (2011a)). Substituting these values 
into Equation (2.14), we obtain 

1.1420 x 10 7 . 



F = 3.7809 x 10 b 



So 



(4.1) 



Starting with the base state described in Section [21 we carry out a modal analysis for 
the turbulent case. The results are shown in Figure [2j Here, we have fixed Re — 2350 and 
e = 0.1, and have chosen the value of m to give a convective/absolute (C/A) transition, 
such that the main saddle point in the contour plot for the largest value of W; has just 
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become positive. The other conditions described in Section 3.1 that are required for 
absolute instability also apply. However, from Figure [2] it is not certain that the system 
is indeed absolutely unstable, given the odd features near the imaginary axis caused 
by different eigenmodes being dominant in different parts of the complex wave number 
plane. 

Therefore, we turn to the ray analysis, representative results from which are shown in 
Figure [3j The two cases are chosen for clarity, since they lie away from (and at opposite 
sides of) the C/A transition such that the respective convective and absolute behaviour 
is clearly visible. In this figure, the amplification of a pulse is very large. This poses a 
practical problem when trying to infer a C/A transition: when starting from an absolutely 
unstable case and repeating the analysis for a slightly different value of the parameter 
m, a convective instability is found when the left tail of the signal decreases as time goes 
by. But the tails of the pulse are close to the part of the ir-domain that is affected by 
numerical error. Moreover, for large to, the strong amplification of all parts of the pulse 
makes a comparison of signal tails difficult. The accuracy with which a C/A transition 
can be determined with the ray analysis is therefore limited. That being as it may be, 
close inspection of the ray analysis results near the C/A transitions inferred from a modal 
analysis confirms the latter. 

Carrying on from these studies, we have constructed a flow-regime map using the 
modal analysis. We have carefully followed the dominant saddle points such as those 
in Figure [2] and have performed the necessary checks required for confirming absolute 
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Figure 3. The norm A(x,t) for (a) m = 400 (convective); (b) m — 1300 (absolute), at t = 30, 
40, 50, 60, 70 (from bottom to topVThe other parameters are r = 1000, Re = 2500, e = 0.1, 
and S and F are given in Equation (4.1 1. The DNS parameters are N\ = 21, N2 ~ 60, L x — 150, 
Ax = 0.01, and At — 0.01. Only a fraction of the x domain is shown. 




Figure 4. Fl ow-r egime map for the turbulent base state. Here, r = 1000 and S and F, are given 
in Equation (4.1|. Mode Ml can be stable (S, neutral curve marked with circles), convectively 
unstable (C), or absolutely unstable (A, C/A transition curve marked with squares). Further 
modes, called M2 and M3, can be convectively unstable, with neutral curves ('NC') given by 
dashed lines and labelled by triangles and diamonds respectively, (a) Variations in e at fixed 
m = 1000; (b) Variations in m at fixed e = 0.1. 



instability (Section 3.1 1. Our modal studies have furthermore been confirmed by the ray 
analysis discussed above. In Figure^a), we see that if Re is increased for a fixed value of 
m, the system generally goes from a stable state, through a convectively unstable to an 
absolutely unstable state. Figure |4^b) shows that the convectively unstable regime can 
disappear entirely, when large values of m are used. The same figure also shows that the 
neutral curve 'pushes' the C/A transition at large values of m back to higher Re- values; 
this is discussed further below. We note in Figure gb) the presence of two additional 
modes of instability whose neutral curves are marked 'NC:M2' and 'NC:M3'; hence, three 
modes are convectively unstable. A standard temporal energy-budget analysis conducted 
at Re = 5000, m = 600, confirms that M2 and M3 both derive most of the destabiliz- 
ing energy from the interfacial region, and a small contribution from the liquid layer. 
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600 




(c) 

Figure 5. The complex wave number along the C/A transition in Figure [4] as a function of 
Re (squares, solid lines). In (a), the wavenumber of the most dangerous temporal mode is also 
shown (circles, dashed line). Filled squares: variation in e at fixed m; empty squares: variation 
in m at fixed e. Panel (c) is a rescaled version of the imaginary wave number in (b). 



The profile of the wave-induced Reynolds stress has also confirmed that they are both 



conventional 'internal modes' in the liquid (Boomkamp & Miesen 1996 O Naraigh et al. 



2011a). However, these modes differ in one respect: M2 has a speed comparable to, but 
less than, the interfacial velocity whereas M3 is much slower, especially at large Reynolds 
numbers. 

In Figure [5] we examine the complex wave number (a r ,a;) at the C/A transition. We 
have plotted a Y as a function of Re for various (e, m) combinations; we have also plotted 
cv^-Reint) 1 / 2 in the same manner. Here Re lnt = r Ui n te/m is the liquid Reynolds number 
based on the base-state velocity at the interface, [/ int (the reason for this rescaling will 
become apparent in what follows). Figure [5^ a, b) demonstrates that a r at the transition is 
governed primarily by Re. The real and imaginary parts of the wave number both increase 
significantly with Re. By rescaling the wavenumbers with respect to the lower-layer depth 
(i.e. letting a r — > ea r ), this variation is given some context: for Re — 1000 the wavelength 
is comparable to the liquid-layer depth, while for Re — 5000 it is approximately 10 times 
smaller than the liquid-layer depth. Figure [5]ja) also indicates that a r at the saddle 
point along the neutral C/A curve closely follows a r for the most-dangerous temporal 
mode. This confirms the observation in Figure[2]that the saddle point lies almost directly 
below the most dangerous temporal mode. This is explained as follows: for all the cases 
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Figure 6. Reduced flow-regime map for the turbulent base state. Filled symbols: variation in 
e at fixed m; empty symbols: variation in m at fixed e. 



considered here, the group velocity duj r /da r calculated for temporal modes depends only 
weakly on a r , and a straightforward application of the Cauchy-Riemann conditions to 
the analytic function uj = u T (a z , a{) + iwi(a r , a{) shows that for such group velocities, the 
a r -values for the most-dangerous temporal mode and the saddle-point mode coincide. 

Most of the results in Figure |4|a,b) collapse into one plot, namely Figure[6ja). Almost 
all corresponding boundaries obtained for Figure |4ja,b) collapse. Especially of interest is 
the large- i?e/large-i?e; nt behaviour of the C/A transition (branch I), which corresponds 
to a critical value of Re m t that is virtually independe nt of Re. It is po s sible t o explain 
this scaling behaviour using the theory developed by O Naraigh et al. (2012) (see also 



Appendix [AJ). There, the spatio-temporal growth rate LJi(a T ,a{) is prescribed in terms 
of a Taylor series in «j, where the coefficients in the Taylor series are derived from the 
purely temporal linear stability analysis, and depend on a x . We have verified that the 
C/A transition curves in Figur e [6| are described accurately by a quadratic truncation 
of this Taylor series (Appendix |A[) . In this 'quadratic approximation', the saddle point 
occurs at a r — a*. , where a* solves 



diO; 



temp ^2^temp 



da r dotf 



-Cg(a*) 



Cfcg 

da r 



(4.2) 



In the current application, the group velocity is approximately constant, hence a* rs 
c*r,max, the location of the most-dangerous temporal mode. The same quadratic approx- 
imation gives the following condition for the onset of absolute instability: 



dy emp 

dai 



temp/ *\ 12/ *\ 



(4.3) 



We now approximate each term in Equation (4.3). Consider first of all the right 



hand side. Previous work ( O Naraigh et aL||2011a | demonstrates that the group velocity 
c g is only slightly in excess of the interfacial velocity U- m t (this is also consistent with 



experiments, e.g. Cohen & Hanratty (19651). Thus, we approximate the group velocity 
as c g ps [7i nt . Moreover, since the liquid layer is thin and viscous, the base-state velocity 
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in the liquid is close to linear shear flow. In addition, for thin films. Re* = 6 Re, w here 9 
is a geometric factor independent of the flow parameters \0 Naraigh et a/.|2011a ). Thus, 
we have the following string of equalities: 

Uvak _ Tjd L _ e_Rel _ ^_g2 Re 
Uq ^lUq m Re m 

But 

PLU int d L r /[/i„t 



RSint — , s , j , - 



hence 

U 



Uint =(e/r l / 2 )Re\> 2 



int ' 



-' ~ {9 2 /r) Re int . (4.4) 



and 

c l 

Next, we consider the left-hand side. We use the quadratic approximation 

Wl tcmp = Aa-|Ba 2 , A,B>0. 

Note that this is not a long-wave approximation, but is instead a fit to the data around 
the most-dangerous temporal mode, where the fitting parameters A and B are selected 
with respect to the actual, computed values of the maximum growth rate, expressed as 
max(o; 1 tomp ) = A 2 /(2B), and the cutoff wave number, written as «o — 2A/B. We now 
consider three parameter regimes and investigate the functional dependence of A and B 
on the Reynolds number Re. 

Functional dependence for fixed e sa 0.1 and varying m: At large Reynolds numbers, 
both the growth rate and the cutoff wave number increase linearly with Re. We write 

max(w 1 temp ) = k\Re, ao — k^Re, 

where k\ and k% are constants of proportionality; these are measured to be m-independent 
(Figure uj). In other words, 

\{A 2 /B) = k t Re, 2A/B = k 2 Re. 

Hence, A 2 = 2k 1 Re(2A/Rek 2 ), or A = Ak 1 /k 2 . At maximum growth, 



temp d"u> ; 



2, temp 



\A 2 = Skl/kl, (4.5) 



2 



da 2 

and the LHS of the C/A transition criterion is independent of Re for large values of Re. 



Combining Equations ( |4.3| , (4.4 1, and (4.5), we have the following criterion for the onset 
of absolute instability: 

lGrkf 



-Rei n t — W, C\ — 



fc 2 6 



where C\ is a parameter-independent constant. Hence, at large Re, there is a critical 
Reynolds number i?ei n t for the onset of absolute instability. This is consistent with figure 
[6] (branch I). 

On the other hand, for smaller values of Re, Figure [7] gives 
max(w 1 temp ) cx \Re — i?e c | p , p > 1, 
where Re c is approximately independent of m. We have also measured the cutoff wave 
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L I'M 



(a) (b) 

Figure 7. Dependence of the temporal problem on the system parameters, (a) The cutoff 
wavenumber ao, for which a; ; emp (ao) = 0; (b) the maximum growth rate max(wj emp ). The 
parameters to be varied are the pair (m, e). Circles: (4000,0.1), Squares: (2000,0.1), Crosses: 
(1000,0.1); Diamonds: (500,0.1); Dashes: (1000,0.005); Stars: (1000,0.05). Trendlines have been 
added to the insets. 



number as 

a = a 0c + k 3 \Re- Re c \ 1/q , q > 1, 

where ao c and k 3 are approximately independent of m. (We have found p s» 5/2, q 
see Figure [7]) Putting these facts together, we get 



temp 



<iV omp 



dai 



1 



n, 



\Re - Re c \ 2p , Re -> Re c 



(Ic 



and the criterion for the onset of absolute instability is therefore 

Re ia t = C 2 \Re-Re c \ 2p , 

where C'2 is independent of m. This relation is consistent with figure [6] (branch II). 

Functional dependence for fixed m and varying e: At large Reynolds numbers and 
relatively large values i?ei nt such that e « 0.1, both the growth rate and the cutoff wave 
number increase linearly with Re (Figure [7|. Thus, the scaling arguments employed for 
fixed e and varying m pertain here also, and the criterion for absolute instability is again 
-Rei n t = Ci, where C\ is a parameter-free constant. This relation is again consistent with 
Figure [6] (branch I). For smaller values of Re, but for still relatively large values of i?ej n t, 
Figure [7] gives 

max(w 1 temp ) cx \Re — i?e c (e)| p , p > 1, 

where i?e c (e) is approximately independent of m but depends on e. We have also measured 
the cutoff wave number as 

ao = a Qc (e) + k 3 \Re- Re c (e)\ 1/q , q>l, 

where aoc(e) is approximately independent of m but dependent on e, and where k 3 is 
independent of m and e. Putting these facts together, we get 



temp 



dV emp 



da 2 



1 



"0e( e ) 



|i?e-i?e c (e)| 2p , Re -> Re c 
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and the criterion for the onset of absolute instability is therefore 

Re int = - C f-\Re-Re c (e)\ 2 P, 

where C\ is independent of m. This relation is consistent with Figure [6] (branch III, large 
-ffei n t), in the sense that the stability boundary depends on -Rej n t, Re, and e. 

As e is reduced, the Reynolds number i?ei n t is also reduced. In such a regime, and for 
large values of Re, we have 

max(w| emp ) = k^eRe 2 jm, ao = k$Re, 
where k^ and k§ are constants of proportionality. In other words, 

\{A 2 /B) = k i eRe 2 /m, 2A/B = k 5 Re, 
Hence, A = Ak^eRe / (k^m) , and B = 8fc 4 e j '{k^m) . Thus, at maximum growth, 



,2, temp 
temp " % 

-W; : 



da? 



2 1.2,™ 2 il - c int- 



'„2 



Re 2 



Therefore, in this case, the LHS of the C/A transition criterion is proportional to Ref nt , 
while the RHS is proportional to Re lnt . Therefore, there is a minimum value of Re m t 
(rather than a maximum value) for absolute instability. This corresponds precisely to the 
lower-branch C/A transition curve in Figure [6] (branch III, small i?e; n t). 

These results have been presented for fixed values of F and So. We briefly sketch 
the effect of varying these parameters. Increasing _Fo is stabilizing, and raises the critical 
Reynolds number Re c for the onset of temporal instability. This suggests that branches 
II and III should shift to the right under such an increase. Similarly, increasing So makes 
the temporal disturbances more stable, the main effect of which is (counter-intuitively) 
to increase — w* cmp d 2 w* cmp /da!r at a T = a rjmax , which implies that the critical Reynolds 
number Re lnt for the onset of absolute instability at large Re and i?ej nt should be raised. 
We have performed some detailed calculations for the onset of absolute instability with 
the 'full' dispersion relation, the results of which agree with this description provided by 
the quadratic approximation. 

We also comment on the scaling behaviour for the imaginary part of the wave number 
at the C/A transition. The value of «i at transition in the quadratic approximation is 

a\ = — 2 max(w 1 ' emp )/c g . 

Approximating c g cx \jRe- m t, we have ai\JRe ln t oc — max(a; 1 tcmp ). For large values of 
the Reynolds numbers Re and i2e m t) we have max(w 1 temp ) cx Re, where the constant 
of proportionality is parameter-free. This gives — ai^/Re ln t oc Re, in agreement with 
Figure [5] 

One further result concerns the alignment of branch I (main C /A transition, Figure [6]) 
with the neutral temporal stability curve of M2. It is as if the absolute instability is 
quenched when M2 becomes unstable. This is due to mode competition, which occurs in 
a convectively unstable regime close to the C/A boundary, in which both Ml and M2 
are unstable, with a wave length comparable to the liquid-layer depth. In Figure [HJ the 
contours of w r are shown for the least stable mode at each complex a. This case involves 
the parameter values Re = 5000, e = 0.1, and m = 600. The spatial curve u>i = of 



the most dangerous mode is identified with Ml by a Gaster-type analysis (see Gaster 



(1962) and Appendix [B|). The contours uj t = Const, that connect orthogonally to the 



spatial curve are identified also with Ml. The orthogonality of contours of the real and 
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(a) m = 600 (b) m = 700 
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(c) m = 710 (d) m = 900 

Figure 8. Spatio-temporal mode competition and the C/A transition: contour plots of u r for 
the most dangerous mode in the complex a-plane. The black lines represent the spatial curve 
Wi = of Ml. Here Re = 5000, and e = 0.1. The M2 temporal mode is unstable for m < 800, and 
the C/A transition occurs near m = 1100. The sharp jumps in oj r represent mode competition 
between Ml and M2, and the saddle point in (a),(b) does not correspond to Ml. 



imaginary parts of an analytic function is a straightforward consequence of the Cauchy- 
Riemann conditions. In this way, we have established in Figure ^a.) that the saddle 
point (a necessary precursor for absolute instability) does not originate from the Ml- 
eigenmode. However, upon increasing m (Figs. J^b-d)), the region of a-space in which 
Ml dominates fills out, and the spatial curve of Ml connects to a saddle: it is as if the 
dominant saddle point 'switches' from non-Mi to Ml. Going to higher values of m, the 
sign of u>i at the saddle point becomes positive, and absolute instability ensues. 

We have investigated in more depth the identity of the rival saddle points by using the 
quadratic approximation, which provides a direct connection between the spatio-temporal 
and temporal modes. By reconstructing Figure [8] from the quadratic approximation (not 
shown), we have confirmed the spatio-temporal mode competition. Care has been taken 
to ensure correct identification of the dispersion curves by double-checking the continuity 
in cv r of the corresponding w r curves. In a reconstruction of the Wi-landscape of Ml, a 
saddle point exists at a r ~ a r ,max- In a further reconstruction, two additional saddle 
points occur, and are connected with M2: one lies near ct rimax , while a further saddle 
point resides at large ce T . Finally, we have verified that as m increases, the values of 
wi associated with Ml increase (especially near the saddle point), eventually leading to 
absolute instability; the other saddles are simultaneously overwhelmed by this increase 
and no longer produce mode competition. 
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4.2. Energy budget 

The purpose of applying an energy-budget analysis is to infer the mechanism for instabil- 
ity generation. In view of the connection between spatio-temporal and temporal modes 
exemplified in the quadratic approximation considered in Section 4.1 a spatio-temporal 
energy-budget analysis necessarily inherits the properties of the corresponding temporal 
analysis. Nevertheless, we consider the energy-budget analysis in this section, and inves- 
tigate how such a study applies to the growth of a pulse (as opposed to a normal mode) . 
The ray analysis provides a means of examining such pulses. 

We examine the mechanism by which the pulse grows, noting that this does not have 
to be the same throughout the pulse. We introduce the kinetic energy density 

Kfat) = \J dzr L \5u\ 2 + \J dzr G \5u\ 2 , (4.6) 

where 5u = (5u, Sw) is the perturbation velocity. We differentiate this expression with 
respect to time at a fixed location x, and apply the equations of motion (2.7 1 and Gauss's 
divergence theorem to obtain the following flux- conservation equation: 

dK dF K , . , . . , , 

-fr+-faT = s B fat) + si fat), (4.7) 

where the source/sink terms sg and s/, and the flux Fk are described in what follows. 
First, we introduce the following notation for the perturbative contribution to the viscous 
stress tensor <5Ty: 

ST ij = -S ij Sp+—^—6u j + —Su, 

we also denote the separate fluid domains by (ol,6l) = (— e, 0) and (og,&g) = (0, do)- 
The flux Fk can then be written as 

F K fat) = f L dz [\r L U (z) \5u\ 2 - 5u5T xx - 6wST xz ] 

J a L 



/ dz [\r G Un (z) \8u\ 2 ~ SuST xx - 5wST xz ] 



In a similar manner, the bulk source/sink term takes the form 
s B fat) = [ REY 3 fa *) + DISSj fa t)} , 



3=L,G 



REYj (x,t) = —rj / dzSudw 

J a; 





dz ' 

b, 



DISS, fat) = / dz 2 (5u x ) 2 + 2 {5w z f + (Su z + 5w x f 

Ke J aj L 

and finally, the interfacial source/sink sj (x,t) is given by 

sj (x, t) = TAN fa t) + NOR fa t) , 
TAN{x, t) = (5T^ X 5u L - 8lf x 5u G ) 



NOR{x, t) = (5T^ Z Sw L - 8Tg 5w G ^ 



1 z=Q ' 
) z=0 ' 

There are no contributions to the energy balance from the perturbation turbulent stresses 
because these are neglected in the quasi-laminar approximation (see Section 2.2). 
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Figure 9. Energy budget for (a,b) Re = 2000, m = 2000 (convective); (c) Re = 4000, m = 300 
(convective); (d-e), Re = 4000, m = 2000. Solid line, TAN; Solid line with circles, REY G ; dotted 
line segments represent absolute value of the same variable as along the rest of the curve. All 
are at t = 50 for r = 1000, e = 0.1. 



As in previous studies of the purely temporal instability, the term TAN is identified 



with the viscosity-contrast instability ('Yih mode') (Yih 1967 Boomkamp & Miesen 
1996). A positive value of this term indicates work done by the perturbations on the 



interface due to the viscosity jump across the interface. Again, in analogy with the purely 
temporal case, the term REYl is due to an instability of Tollmien-Schlichting type in 
the bulk liquid flow, while positive values of REYq correspond to an instability of the 
Miles type (Miles 1957) near the critical layer in the gas. Equally, the terms REY L and 
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REYq can be thought of as giving the rate at which energy is transferred from the mean 
flow to the disturbance via the wave-induced Reynolds stresses. 

We examine these terms in detail now (the sinks DISSl.g ar e not of interest, since 
they are necessarily stabilising). To this end, we have selected three states (i?e,m) from 
Figure Qb) that intersect the near-horizontal and near- vertical parts of the C/A tran- 
sition. For these states, we plot the sources and sinks as a function of x at a fixed 
point in time (the t-value is chosen such that all transience has been eliminated from 
the pulse). Naturally, the curves exhibit oscillatory behaviour (the distribution of the 
phase of even a single temporal mode shows oscillations). The energy budget for the 
purely temporal study is averaged over a single wave length, but spatial averaging is not 
conducted here, since the spatial distribution is the focus of the study. Therefore, we 
examine in Figure ^ a) a snapshot of the spatial distribution of the largest terms in the 
budget for Re = 2000, m = 2000, e = 0.1. This parameter set is seen in Figure [4](b) to 
be convectively unstable and to lie to the left of the C/A critical curve. The term TAN 
dominates throughout the pulse, followed by REYq (other, smaller terms are not shown). 
For Re = 4000, m = 300, e = 0.1, this is also the case but, although TAN reaches locally 
the largest values over a wave length, its sign changes, whereas REYq is positive virtually 
throughout the pulse (a zoomed view is shown in Figure gc)). Finally, crossing the C/A 
boundary (Figure [9jd-e) ) causes the distribution of REYq to become asymmetric, while 
the features regarding the signs of TAN and REYq previously observed for m = 300 
still hold. 

Having characterized the turbulent base state in detail, we revisit the laminar problem 
and investigate whether large-r absolute instability is possible there. 



5. Revisiting the laminar problem 

In this section we review the problem of interfacial instability where the upper layer 
is laminar. Although the main conclusions for the turbulent case carry over here, some 
differences arise. The base-state flow Uo(z) is determined in a standard fashion by solving 
the momentum balance 

d 2 U Q _ dP 
dz 2 dL 



and is subject to continuity of velocity and shear stress across the interface at z = 0. To 
facilitate comparison with previous work on laminar flows, for this section only, we adopt 



the non-dimensionalization scheme of Valluri et al. (2010). We set e = dj^jidh + da) and 



Re = pqV (dL+dc) I hg, where the characteristic velocity V is chosen to be the superficial 
velocity (d^ + dc)~ 1 Uo(z)dz. The gravity and surface tension are parameterized as 
Q '■= (pl — PG)g{dL + da) 2 /(pqV) and S :— j/(ij,gV) respectively; the value of these 
parameters is varied in the following parameter studies. 

Because the model for the turbulent case relies on the quasi-laminar theory for the Orr- 
Sommerfeld perturbation equations, the streamfunction equations and the transient-DNS 
numerical method carry over directly to the laminar case, once the necessary rescaling 
and non-dimensionalization have been performed (e.g., using Q and S in the normal stress 
condition) . 

5.1. Parametric study (1): Q and S taken as constant 



We set Q — and take S = 0.01, corresponding exactly to the paper of Valluri et al. 



(2010); this choice enables us to relate the current investigation to the results of Valluri 



et al. (20101 for liquid/liquid systems. For these parameter values, the ray analysis has 
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(e) e = 0.18 

Figure 10. Snapshots of the norm A(x, t) for different values the depth of e for the laminar 
base state. The parameters are Re = 1500, m = 150, r = 1000, S = 0.01, and G = 0. (a) 
t = 100,200,300,500,700; (b)-(d) t = 100,200,300,400,500; (e) t = 50,100,150,200,250. Fig- 
ures (a) and (e) are convectively unstable cases; (b)-(d) are absolutely unstable. In all cases, 
L x — SO, Ni = 21, A^2 =51, the timestep is At = 0.1, and the grid size Ax in the x-direction is 
0.003. 



revealed that the laminar base state is absolutely unstable for a significant portion of 
parameter space, even for a density ratio of r = 1000. This is seen in Figure 10 where the 
L 2 -norm A(x,t) associated with the broadband disturbance is shown at several times. 
These curves have been truncated at the points where they decrease below a specified 
fraction of their maximum value, as the finite working precision makes the results too 
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Figure 11. Flow-regime map showing laminar base state in parametric study (1). Here, 
r = 1000, S = 0.01, and Q — 0. The system is always either convectively (C) unstable or 
absolutely (A) unstable. A second temporal mode is also unstable below the neutral curve 
(dashed curve labelled 'M2', with triangles), (a) Variations in e at fixed m = 150; (b) Variations 
in m at fixed thickness e = 0.1. 
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Figure 12. Flow-regime map of parametric study (1) in terms of the liquid Reynolds number 
i?ei n t. Filled symbols: variations in e at constant m = 150; open symbols: variations in m at fixed 
e = 0.1. As described in the text, a further mode (convectively unstable) exists for Re > 8000, 
whose neutral curve is unaffected by changes in -Reint, m, and e. 



uncertain far away from the pulse. The existence of absolute instability at large density 
ratios was not found in the modal analysis of Valluri et al. (20101, who only ascertained 
that density-matched fluids can become absolutely unstable. The cause for this oversight 
is that the magnitude of the wave number at the saddle point is very large, beyond the 
range searched by Valluri et al. (20101. 

Motivated by the plots in Figure 10 showing absolute instability, we revisit the other- 
wise more accurate modal analysis and perform a large scan through the complex a-planc 
to obtain the C/A boundaries. These boundaries are consistent with the results of the 
ray analysis shown in Figure [TUj Also, the results are similar to those in Valluri et al.\ 
(2010) for r = 1, suggesting that the transition found by Valluri et al. (20101 extends to 
density ratios of at least r — 1000. We have verified that, upon lowering the density ratio 
at a point in the absolute regime of Figure [HJ also identified as such in Valluri et al. 
(2010), the system remains absolutely unstable at intermediate values of r. In contrast 
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Figure 13. The real (a) and imaginary (b) component of the wav e number at the saddle point 
along all C/A transitions in Figure 11 a) (squares) and Figure [TT^ b) (triangles), as functions of 
i?ei n t. The open symbols represent secondary saddle points that do not correspond to absolute 
instability. 



to the r = 1 case in Valluri et al. (2010), for r = 1000, absolute instability occurs only 
at large m- values. 

As in the turbulent case, we have examined the flow-regime boundaries in the (Re, Rei n t) 
plane, where Re- m t :— reReU- m t/m is a Reynolds number based on the liquid-film proper- 
ties and the interfacial velocity U- irA of the base state. In contrast to the turbulent case, 
the results of Figure [FL] do not collapse, although the overall trends otherwise bear some 
resemblance to the turbulent case (e.g. Figur e |4| : b ranch III corresponds to a critical 
value of i?e; n t (consistent with Valluri et al. (2010) for r = 1), while branch II corre- 
sponds to a critical value of Re. These two branches shift when changing the value of m, 
respectively e, thereby contracting or expanding the absolutely unstable regime. Instead, 
we plot the regime boundaries in the (Re, Re- ln t/e) plane in Figure 12 Although the two 
branches labelled 'IF virtually coincide in this plane (this is discussed further below), 
the functional form i?ei n t = ef(Re) of the neutral curve differs substantially from that 
already encountered in the turbulent study. 

We also investigate variations in the wave number at the saddle point (taken at the 
onset of absolute instability) in Figure 13 where e(a T ,a{) is plotted as a function of 
i?ei nt . The contrast between these results and the earlier turbulent results is remarkable. 
The wave number is 0(l/e). While the flow-regime map does not collapse in (Re, Rei nt )- 
space, the wave numbers e(cv r , a{) depend mainly on Re m t, and not separately on Re. The 
wave number «o at the saddle point follows from (duj/da) ao — 0, which therefore only 
involves a relatively simple scaling, whereas the flow regime is determined by Wi(ao) > 0, 
which involves more parameters. Nevertheless, the dependency of ea on i?ei n t is complex. 
Above Re mt ~ 5 (corresponding to branch I) the curves nearly collapse; a small differ- 
ence remains, as is also seen in the corresponding conditions C/A transition curves in 
Figure 12 Upon decreasing the value of i?ej n t, two saddle points emerge. However, only 
one of the saddle points produces absolute instability. Most of the results in Figure |l"3| lie 
on the branch that corresponds to somewhat longer waves (note, the data represented by 
the filled triangles in Figure [13^ a) jump from the upper to the lower branch when i?ei nt 
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Figure 14. Flow-regime map for parametric study (2). Filled symbols: variations in e at 
constant m = 100; open symbols: variations in m at fixed e = 0.1. 



is decreased). These substantial differences between the turbulent and the laminar cases 
are now addressed in a further parametric study. 

5.2. Parametric study (2): Q and S vary inversely with Reynolds number 

We extend parametric study (1) to allow for effects of gravity and surface tension cor- 
responding to systems with a large density ratio, and we set Q — Q (r — I) /Re and 
S = So/ Re. This extension is more reflective of real systems because it corresponds to 
a fixed value of the dimensional surface tension 7, while allowing for variations in S 
through changes in the Reynolds number. Here, the values of Go and So are the same as 
F a and So (Section |4| but such that these are for channels that are ten times smaller: 
£0 = F /10 3 and So = So/10 unless indicated otherwise (the effect of changing these 
values is discussed below) . The results of this parameter modification are shown in Fig- 
ure 14 A striking feature of this study is the blurring of the previously sharp distinction 



between the turbulent case (Figure rc5j), and the laminar study (Figure 12). Moreover, 
increasing Q and S leads to a substantial modification of the absolute region in param- 
eter space. The role played by Q and S in this modification is further highlighted in 
Tables [3^(4] Tabic [3] corresponds to branch I, and demonstrates that an increase in Q or 
S at fixed m and e calls for an increase in Re to sustain absolute instability. Similarly, 
Table [3] corresponds to branch II, and demonstrates that an increase in Q at fixed Re 
and e requires a corresponding, destabilizing increase in m to sustain absolute instability. 
This produces a proportional decrease in Re int . The situation concerning increases in S 
is somewhat counter-intuitive (increases in S produce decreases in the critical value of 
m, which produce increases in the critical value of i?e; n t). However, this agrees with the 
turbulent case, which has already been explained using the quadratic approximation. 

The scaling of the transition curves in Figure [14] can again be elucidated with the 
quadratic approximation. For example, for large Re and i?e; n t, we have measured 

max( Wl tomp ) cx (e/m 1 / 2 )^ 1 / 4 , a - a 0c oc (e 1 / 2 /m 1 / 2 )i?e 3 / 4 , c g « U int oc e/m. 

(s.i) 

Plugging these scaling rules into the quadratic approximation gives eRe/m 2 = const., 
or i?ej n t/e = const, at large Re and i?ei n t, which implies the existence of a critical 
value of i?ei n t /e along branch I. The dramatic distinction between these scaling rules 
and those encountered in the turbulent case owes not to some deep distinction between 
the respective base states, but rather arises from the distinct non-dimensionalization 
schemes chosen in each case. Indeed, as a final test, we have verified that a laminar 
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Qo/Qo,vel So/S(), rcf R e -Rein 



0.1 
1 
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1 
1 
1 



1 
1 
1 

0.1 

1 

100 



2050 
3450 
9250 
2750 
3450 
5350 



0.17 
0.28 
0.76 
0.23 
0.28 
0.76 



1.8-0.65i 
2.2 - 0.9i 
3.2 - 1.81 
4.1 - 2i 
2.2 - 0.9i 
1.2 - 0.45i 



Table 3. Dependence of the critical parameters on surface tension and gravity at low iteint 
(branch I, m = 1000, e — 0.1). The wavenumber at the saddle point is also stated. The un- 
certainty in Re is ±50, which also introduces an uncertainty in the values of Re^t and a. The 
symbols (Jo.rcf and 5o, r of refer to the reference values discussed at the beginning of the section. 
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Table 4. Dependence of the critical parameters on surface tension and gravity at low Re (branch 
II, Re — 10 4 , e = 0.1). The wavenumber at the saddle point is also stated. The uncertainty in 
m is ±2.5, which also introduces an uncertainty in the values of Reint and a. 
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Figure 15. The real (a) and imagi nar y (b) component of the wav e number at the saddle point 
along all C/A transitions in Figure 11 a) (squares) and Figure [TT| b) (triangles), as functions of 
Re ln t- The open symbols represent secondary saddle points that do not correspond to absolute 
instability. 



temporal study, based on the non-dimensionalization scheme previously applied in the 



turbulent case, yields the scaling rules described in Section 4.1 (not shown). 

In Figure 15 we examine again the wave number e(a r ,a;) at the saddle point, at 
the inception of absolute instability. The contrast between this figure and parametric 



study (1) is remarkable. In Figure 15 the results resemble those for turbulence, with one 
significant difference: the saddle-point mode corresponds to a much longer wave. The 
reason for the contrast between parametric studies (1) and (2) is due to the fact that 
the wave length is controlled mostly by surface tension (as demonstrated by Table [4]) ; 
the surface tension in study (2) greatly exceeds that in study (1), consequently, the most 
unstable wave length is longer. 

Referring to parametric study (1) and figure [TT^b) , branch II (varying m) of the C/A 
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Figure 16. Competition between spatio-temporal modes for the laminar case. Shown are con- 
tour plots of u T for the most dangerous mode in the complex a plane at Re — 20000, e = 0.1 for 
different values of m. The black lines represent the spatial curve Wi = of Ml. The sharp jumps 
in uj t represent mode competition between Ml and M2. The M2 temporal mode is unstable for 
m < 75, and the C/A transition occurs near m = 170. 



transition 'follows' closely the neutral stability curve for a second temporal mode, M2. 
The same neutral curve is also present in parametric study (2), but is further removed 
from the C/A transition curve, and occurs at Re- m t/e s» 1.5 x 10 4 . As in the turbulent 
case, it is as if the existence of the M2 mode quenches the Ml absolute instability. Again, 
as in the turbulent case, there is clear evidence of competition between spatio-temporal 



modes, in addition to temporal mode competition. In Figure 16 cj t is shown for the 
least stable mode at each complex a. Results are presented at different values of m for 
Re — 20000, e = 0.1, and r = 1000. The value m = 85 describes a situation where M2 has 
just become stable while Ml is convectively unstable (the point of neutral stability for M2 
is at m = 75). The C/A transition for Ml occurs at m = 170. The spatial curve u>i = of 
the most dangerous mode is identified with Ml by a Gaster-type analysis ( Gaster|[l962 ). 
The contours w r = Const, that connect orthogonally to the spatial curve are identified 
also with Ml. In this way, we have established that the saddle point in figure [l6|[a) 
does not correspond to Ml. Thus mode competition interferes most dramatically with 
the saddle point near the point of neutral temporal stability of M2 (near to = 75), 
while an approach to the C/A transition from within the convective regime (i.e., an 
increase in to from a low value) causes the mode competition to disappear gradually, 
thus producing a conventional single-mode saddle point. This strongly suggests that the 
proximity of the M2 neutral curve and the Ml C/A transition in the flow diagram is no 
coincidence, and is in agreement with the earlier results for the turbulent case. The same 
phenomenon persists for lower Reynolds numbers (e.g. Re = 5000). However, we have 
presented the results for Re = 20000 because this regime exhibits the mode competition 
most clearly. At such high Reynolds numbers, a third unstable mode comes into existence 
(visible in Fig. 16 in the neighbourhood = and a r = 2). This mode has a critical 
Reynolds number Re rs 8000 and its neutral curve is almost independent of i?ei n t, to, 
and e. However, this third mode does not play any role in the spatio-temporal mode 
competition and is not discussed any further. 

Finally, we have applied the spatio-temporal energy budget to the laminar flow for the 
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e <w(ST) A max (Modal) a(ST) Q(Modal) 



0.18 


0.0934 


0.0936 


15.2 


16.3 


0.15 


0.0865 


0.0866 


16.9 


16.5 


0.12 


0.0787 


0.0788 


17.6 


17.2 


0.08 


0.0633 


0.0633 




20.8 


0.04 


0.0326 


0.0326 




30.9 


0.025 


0.0179 


0.0179 




43.3 



Table 5. Comparison between modal and spatio-temporal analyses at various values of e, at 
r = 1000, m = 150, and Re — 1500. The rest of the parameters are the same as in Figure [l0"| 



cases corresponding to Figure 10 'a), (c) and (e). The results indicate that over most of 
the pulse width, TAN dominates: this represents a source of energy for the instability 
that is associated with the viscosity-contrast mechanism, which dominates in temporal 
instability. However, REYq also plays a significant role. In all three cases, REYq is 
negative on the downstream side of the pulse, but positive on the upstream side. This 
is most visible in the absolutely unstable case. In contrast, for the turbulent base state, 
it is only in the absolutely unstable case that REYq develops this asymmetry, where it 
even becomes the dominant term on the upstream side in the energy-budget analysis. 



5.3. Ray analysis revisited 
Some further information about the absolute instability of Ml is obtained from a ray 



analysis (Figure 17). The parameters are taken from the cases studied in Figure 10 The 
branches of the C/A boundary in Figure [TTJa) can be mapped on to families of o~(v)- 
curves in Figure 17 decreasing e from e « 0.18 in the latter gives a family of cr(z;)-curves 



that is associated with the upper branch I of the C/A transition in Figure [TT[a) (m 
constant, e varying). As e is decreased, the £r(u)-curves shift downwards in their entirety, 
until the growth rate o~(v = 0) eventually becomes positive, indicating a switchover to 
absolute instability. Conversely, increasing e from e 0.025 gives a family of cr(w)-curves 
that is associated with the lower branch III in Figure [TT^a) . These curves are steep for 
small e-values, such that o~{v) is negative. As e increases, the slope diminishes, until 
a(v — 0) is positive, and absolute instability is attained from below. Furthermore, in a 
convectively unstable regime at large e-values, cr(v) is linear in v for small v, and the 
corresponding a T - and Wi-values approximate negative constants, consistent with Equa- 
tions (3.4). Again with reference to Equations (3.4), the upper branch in Figure 11 'a) 
corresponds to a shift in the value v for which a = cr max (hence ot\ ~ 0), whereas the 
transition at the lower branch is associated with an increase in temporal relative to spa- 
tial growth. Finally, for large values of e, the cr(ti)-curve is very wide, indicating a rapid 
spreading of an initial pulse, whereas for low values the pulse remains more localized. 

For the cases considered herein, we have verified that the wavenumber of the spatio- 
tcmporally most dangerous mode coincides with that of the temporally most dangerous 
mode (Table [5|. It was not possible to compute this wavenumber for all of the e-values, 
and our finite-difference discretization of Equation (3.4d) is not always accurate. It is 



for this reason that we did not include a r -v plots in Figure [TTj Nevertheless, Figure [17] 
still contains important information: a previous set of purely spatial and purely temporal 
analyses (see Valluri et al. (2010) and Appendix [B]) found extremely large purely spatial 
growth rates (exceeding the least-negative purely spatial growth rate by an order of 
magnitude) that were not excited in a DNS. From the spatial growth rates along rays 
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Figure 17. Ray analysis: (a) spatio-temporal growth rates, (b) spatial growth rates along rays; 
(c) temporal growth rates along rays, for laminar flow at various values of e, at r = 1000, 
m = 150, Re — 1500. The rest of the parameters are the same as in Figure [To] Trend lines have 
been added in (a) merely to guide the eye. 



in Figure [l7|b) , we see that such extremely large spatial growth rates are not selected, 
instead, the spatial growth rate comes from the same eigenmode as the temporally most 
dangerous mode. 



6. Discussion 

6.1. Flow-regime maps 

It is of interest to compare our regime boundaries (e.g., Figure |6| with flow-regime maps 
determined from experiments. Our theory is of course applicable to water/air systems, 
for which channel-flow experimental studies are available. However, the results in Sec. [4] 
prove that water/air systems are usually convectively unstable, while the main focus 
of the present paper is on absolutely unstable systems. Consequently, we have searched 
the literature for experiments concerning air and oil, for which the viscosity contrast is 
much higher, and for which absolute instability is expected. We have found that such 



experiments were conducted on pipes (with the exception of the work by Gondret et al. 



(1999) on Hele-Shaw cells). Nevertheless, a comparison between our channel-flow model 
and a pipe-flow experiment is justified, since a 2D analysis may be representative of the 
mid-plane in pipe flow, especially for thin films, linear transverse waves (in the third 
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Figure 18. Flow-regime map Figure [4] in terms of superficial gas and liquid velocities. The 
curves with symbols are the present work; other solid lines and classifications in italics are 
experimental findings by Hoogendoorn (1959) for spindle oil/air flow in a 91mm diameter hori- 
zontal pipe wherein the following tiow regimes are identified: stratified flow with a fiat interface 
(S), wavy-stratified flow (W-S), wavy-stratified flow with suspended droplets (MIST), and slug 
flow (SL). 



dimension) are usually overwhelmed by 2D waves, and in any case it is necessary to 
demonstrate that the boundaries determined in our study are not restricted to an obscure 
limit. 

The results of Figure Qa) are therefore presented in Figure 18 in terms of superficial 
velocities Uq and Ul, where c g = Ug/{Ug + Ul) and U m = Ug + Ul is the mixture 
velocity for parameter values approximately corresponding to the flow of air and spindle 
oil through a pipe of 91mm diameter, and flow-regime boundaries determined experimen- 
tally by Hoogendoorn ( 1959[ ) are superposed (our data are for m = fOOO, r = fOOO and 
the choice of S and F mentioned at the start of Secti on [4} w hich are also representative of 
the experimental conditions stated in Hoogendoorn (1959)). The more widely used flow- 
regime map by Mandane et al. ( 1974 1 cannot be used here, as it is for water/air systems, 



but they observed good agreement with water/air data also published by Hoogendoorn 
dl959| . 

We first note from Figure[f8]that the neutral stability curve from the temporal analysis 
agrees reasonably well with the experimental data on the transition from stratified flow 
(i.e., no waves) to wavy stratified flow. This is expected b ased on the comparisons with 
experiments in our previous work on temporal instability, O Naraigh et al. (2011a I. We 



note here that the superficial velocities in a 2D system are expected to be larger than 
for a pipe flow at the same maximal velocities in a fluid, which may have reduced the 
difference between the results somewhat. On the other hand, the fact that the liquid 
layer at the centre of the pipe is expected to be smaller due to gas pushing the liquid 
to the side of the pipe is expected to be of less significance, as the liquid height is the 
parameter along the theoretical data in the figure. 

The main C/A transition is predicted theoretically here at a total superficial velocity 
somewhat below the experimentally determined transition from wavy-stratified flow to 
mist flow (waves with droplets). A relation with the onset of atomization for very viscous 
liquids cannot be ruled out. In fact the experimental mist flow regime lies almost entirely 
in the absolutely unstable regime, since the critical value of the parameter c g for this 
regime is close to the high-i?e C/A transition from the theory (and even more so, the M2 
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neutral curve). In the experiments, slug flow is observed at c g below this critical value at 
large U m . Hence slug flow is observed experimentally when according to the 2D theory 
the flow is convectively unstable. 

6.2. Potential use in global and non-linear analyses 

Although we have identified and classified the nature of the linear instability for our 
stated two-phase stratified problem, we briefly outline extensions that follow naturally 
from the present work. 

The stability analysis considered throughout this paper is a 'local analysis' in the 
sense that the base flow varies on a length scale which is long relative to the wave 
length of the instability waves (this being true whether the analysis is linear or not). 
A more general approach is global stability analysis, wherein the whole of the physical 
domain is considered (Huerre & Monkewitz 1990 Chomaz 2005). Under such a technique, 



spatial streamwise variation is accounted for in both the base flow and the perturbation 
terms, permitting the study of nonparallel flows, i.e., no restrictions are placed on spatial 
scalings. Of course, such studies require appropriate computational power and therefore 
have only been pursued recently. In the context of the present work, such an approach 
permits changes in the thickness of the liquid layer which is important for many of the 
stated applications, but it means that ray analysis in the form used herein can no longer 
be employed, as periodic boundary conditions no longer apply. 

The theoretical analysis required to elucidate non-linear effects discussed herein is 
elaborate. Thus, an important additional step might involve a numerical (e.g., DNS) 
study and experiments, with the aim of identifying which key non-linear effects dom- 
inate. Specifically, whether a transition to absolute instability in the nonlinear regime 



may precede the linear C/A transition (this has been observed by Gondret et al. (1999) 
in Hele-Shaw cells). Another question that a full numerical simulation might address 
concerns the location of the dominant non-linear interactions, i.e., at the phase interface 
or in the bulk of the liquid layer. A full numerical simulation would also permit a com- 
parison between periodic boundary conditions (as used in the present work), and inlet 
conditions (as in a true channel flow). In this paper, the ray analysis is conducted in 
Fourier space, and required periodic boundary conditions. However, the pulse-type ini- 
tial conditions considered herein involve extremely long domains, and comparisons with 
even larger domains are self-consistent. Nevertheless, full numerical simulations would 
further validate these observations. Such a study and the subsequent non-linear analysis 
will form the basis of future work. 



6.3. Summary 

We have studied the linear stability of a liquid layer sheared by laminar or turbulent 
gas flow in a two-dimensional Cartesian geometry. A Poiseuille (laminar) or Reynolds- 
averaged (turbulent) model is used to describe the stratified two-phase base flow. Results 
from a normal-mode analysis and a ray analysis show that the base flow is absolutely 
unstable to linear perturbations for large regions of parameter space. In particular, for 
density ratios of 0(1000), clear evidence is given of oil/gas flows being absolutely unstable 
for a sufficiently large dynamic viscosity ratio (at least O(10 2 ) for the laminar and O(10 3 ) 
for the turbulent base state), provided the Reynolds number of the flow is sufficiently 
high. Since the dynamic viscosity ratio for air/water systems is only m ~ 55 at 20°C, the 
present results for the turbulent case lead us to conclude that laboratory experiments 
carried out for water/air may not be representative for oil/gas systems. 

In both turbulent and laminar base states, the flow-regime map collapses in the (Re, 
Rei n t) plane, where Re lnt is the liquid-based Reynolds number (with a small subset of 
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exceptional cases). A recently developed theory for C/A transitions (O Naraigh et al. 
|2012[ ) has enabled us to formulate a simple criterion for the onset of absolute instability, 
based on a competition between instability growth and convection by the group velocity. 
A single number, based on the temporal growth properties, encodes the tendency for 
instability growth. This number is governed by e, m, and Re, but becomes parameter- 
free in certain limits, while the behaviour of the group velocity is governed simply by 
i?ej nt . This criterion therefore explains the importance of Re mt in the C/A transition 
curves, and, moreover, when examined in detail, explains the collapse of the transition 
curves into universal forms. 

Independent confirmation of our results from the modal analysis has been obtained us- 
ing a ray analysis. The ray analysis benefits from our development of an efficient method 
of solution for linear differential-algebraic equations, which allows us to extend dras- 
tically the simulation time over which the evolution of an initially pulsed disturbance 



can be traced, thereby eliminating the problems encountered with DNS by |Valluri et al. 



(2010). This method all but guarantees the detection of absolute instability, if it exists. 
Additionally, the ray analysis leads naturally to an energy-budget method that deter- 
mines the origin of the spatio-temporal instability. Results from this method demonstrate 
the importance of the viscosity-contrast mechanism acting at the interface, as well as a 
transfer of energy from the bulk gas flow to the waves, as the sources of the instability. 
Nevertheless, there exist several internal modes (which also derive energy from the in- 
terface) that are at least as significant in other cases. A further application of the ray 
analysis is envisaged for flows wherein the complex dispersion relation contains singu- 



larities (e.g. if ui(a) has a pole (Healey 2007 20091 or root- type behaviour (Aships & 
Reshotko"fT9 90 ) at some point ao). In such problems, the application of the saddle-point 



criterion, with its associated construction of the steepest-descent path enclosing all the 



dynamically-relevant singularities, is a non-obvious task ( |Rees fc Juniper 2010[ ). For these 



cases, a straightforward ray analysis could be used to verify the correctness of one's im- 
plementation of steepest-descent method. The ray analysis suffers from the amplification 
of numerical error when the spatial growth rates are large, meaning that only a small 
region around the impulse centre can be used to extract meaningful information. Thus, 
both methods are necessary to obtain a complete picture of the stability properties. The 
ray analysis and the quadratic approximation are quite general, and we anticipate that 
a wide range of multiphase flow scenarios can be tackled using these techniques. 



Appendix A. Validity of the quadratic approximation as applied to 
the C/A transition of two-phase stratified flow 



We review the so-called 'quadratic approximation' developed by O Naraigh et al. 



(2012) and demonstrate that it approximates well the precise criteria for the onset of 
absolute instability. The quadratic approximation is based on the following identity for 
the analytic continuation of the growth rate into the complex plane: 



(-i)" d^ Cg 2n+1 , ~ d 2 »+y omp 



I \ temp/ \ . \ 1 " L g 2n+l i 



^(2n + l)\da 2n 1 ^(2n + 2)\ da 2 , 

n=0 n—Q x 7 



2n+2 



a 



2n+2 



(Al) 

where c g is the group velocity dw r /da r in a purely temporal analysis. This is a conse- 
quence of the Cauchy-Riemann conditions on uj{a) viewed as a holomorphic function 
on an appropriate open su bset in the complex plane, and has been derived elsewhere 



by O Naraigh et al. (2012 1. We truncate this series at quadratic order in a>\ and apply 
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(a) (b) 

Figure 19. Comparison of flow-regime maps resulting from the quadratic approximation (open 
symbols, keeping m constant; filled symbols, keeping e constant) and the full modal analysis 
(solid lines) for the turbulent (a) and laminar (b) base state. 



the conditions du>i/da Y = daJi/dcti = for a saddle point. This yields the conditions 

Simultaneous solution of these equations yields the following a r -value for the location of 
the saddle point: 

c(a) dc B _ rf^ cmp rfy° mp | i c 2 (a ) ^ 3 ^ tcmp / rfV cmp \ (A3) 

s r da r da r do^ 2 s r \ daf / da 2 I 

(Note that for a strictly quadratic approximation, as in the main part of the paper, the 
third derivative should be neglected here.) At the onset of absolute instability, = 
at the saddle point. Substitution of this condition into the quadratic approximation to 



Equation ( A 1 ) yields the following criterion for the onset of absolute instability: 

,2 temp 

5 temp / * \ 1 2/ *\ m a 

w ; p (a r ) = 2C g (a r ). (A 4) 



doti 



where a* is the root of Equation (A3) 



In Figure 19 the results of the quadratic approximation are compared with the full 
modal results that have been taken from Figures [6] and [14] for the turbulent and laminar 
base states, respectively. All branches are qualitatively correctly predicted, and quanti- 
tatively accurate results are even obtained at some, especially all results at large Re in 
the turbulent case, and those for low Re and a large viscosity ratio in both laminar and 
turbulent cases. 



Appendix B. Analysis of purely spatial modes 

The focus throughout the paper has been on spatio-temporal analysis. In particular, 
our modal analysis has involved the extraction of saddle points from the solution of 
an eigenvalue problem in a complex wave number space. In this section, we review the 
Gaster theory for the extraction of purely spatial growth rates. The focus herein is on 
the turbulent base state. 

The most basic description of linear instability is a temporal analysis, which involves 
the solution of the eigenvalue problem (2.15) for a = a r only. This gives a dispersion 
relation (w r (a r ), Wi(a r )). The pair (w r (a r ), Wi(a r )) that maximizes Wi is called the most 
dangerous temporal mode. The maximum is taken over the complete spectrum of modes 
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a a 
(a) (b) 

Figure 20. (a) Spatial stability study; (b) Comparison with Gaster analysis. Here 'MDF' 
refers to the 'most-dangerous frequency', or, equivalently, the most-dangerous mode. 



in the linear Orr-Sommcrfcld (OS) equations. The system is unstable if cjj > at the 
most dangerous temporal mode. This case has been investigated e xhaustively by the 
present authors ( |6 Naraigh & Spelt|[201o| |Q Naraigh et aZ.||2011a.|6 ) . 

We therefore turn to the spatial case, wherein we are concerned with the dispersion 
relation ai = a x (a r ,Wi = 0), which can be obtained by treating the OS problem (2.151 



as a non-linear eigenvalue problem in a. In practice, we do not solve this complicated 
problem. Instead, we solve Equation ( 2.15[ ) in the complex a-plane, a :— a r + ictj. The 
result is a discrete spectrum of frequencies {w n (a)}„, which we order according to the 
magnitude of 9 [u„(a)]: the most dangerous mode at complex wave number a is denoted 
by ujc(oi), and is defined such that 



^s[u3c{a)]= max \ 9? [u n (a)} \ 

n. all modes L it 



We plot the zero contour of ^s\uj c {oi)\ = 0, which gives the desired relation = 
a; (a r ,Wi = 0). We repeat these steps for the less dangerous modes. This results in multi- 
ple spatial curves. These are displayed in Figure [20^ a) for the parameter values m = 55, 
r = 1000, e — 0.05, with F and S given by Equation (4.1 ). However, in this figure, we plot 



only the spatial curves produced by the most dangerous mode; the less dangerous modes 
do not produce the same large spatial amplification (spatial amplification corresponds to 
a?i < 0), and are not shown. The large spatial growth rates in the figure would appear 
to dominate in an evolving flow. In practice, however, they are not observed, neither in 
turbulence nor in laminar flows. 

Because the large spatial growth rates are not accessed, we focus our attention on the 
spatial curves near ci{ — 0, for which an analytical theory is available. Gaster ( 1962 ) 
has shown that small spatial growth rates are related to the temporal ones through the 
formula 

/dui T 

- «i = Wi / - — , 



(Bl) 



where the quantities on the right-hand side are derived from the purely temporal analysis. 



We test the applicability of Equation ( B 1 1 in Figure 20 b) . We plot the Gaster curves 



associated with the temporally most dangerous mode and the second most dangerous 
mode in the figure, and compare the result with the contour in Figure 20 '&). There is 



good agreement between the two theories. However, the curve generated from the contour 
analysis exhibits 'kinks' where the most dangerous mode switches from one eigenmodc 
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to another. To obtain perfect overlap between the two curve sets, it is necessary to pick 
out all the eigenmodes from the contour analysis and to cut and re-connect the contours 
so as to obtain smooth curves. This exercise is difficult but yields no new information, 
and we do not pursue it here. Instead, we focus on the Gaster study, where two spatial 
curves overlap with the contour-generated ones. The Gaster curves are obtained from 
the first and second temporally most dangerous modes. Thus, the character of a spatial 
instability is identified with the character of the temporal instability through the Gaster 



formula ( B 1 1 . The Gaster curve suggests that the second most dangerous temporal mode 
contributes most strongly to the spatial instability. However, it is not obvious from this 
simple analysis which eigenmode is excited in an impulse-response scenario: this question 
is studied in detail in Section[5] where it is demonstrated that the spatial curve associated 
with the temporally most dangerous mode is the one that is selected. 

Acknowledgements 

The authors would also like to thank J.-C. Loiseau for carrying out preliminary numerical 
investigations for this project. 



REFERENCES 

Ar, S. & Cai, J.-Y. 1994 Reliable benchmarks using numerical instability. In Proceedings of 

the fifth annual ACM-SIAM Symposium on Discrete Algorithms . 
Arfken, G. B. & Weber, H. J. 2001 Mathematical methods for physicists, 5th edn. Harcourt. 
Aships, D. E. & Reshotko, E. 1990 The vibrating ribbon problem revisited. J. Fluid Mech. 

213, 531. 

Biberg, D. 2007 A mathematical model for two-phase stratified turbulent duct flow. Multiph. 
ScLTech. 19, 1. 

Boomkamp, P. A. M., Boersma, B. J., Miesen, R. H. M. & v. Beijnon, G. 1997 A Cheby- 
shev collocation method for solving two-phase flow stability problems. J. Comp. Phys 132, 
191. 

Boomkamp, P. A. M. & Miesen, R. H. M. 1996 Classification of instabilities in parallel 
two-phase flow. Int. J. Multiphase Flow 22, 67. 

Bradshaw, P. 1974 Possible origin of Prandtl's mixing-length theory. Nature 249, 135. 

Chomaz, J.-M. 2005 Global instabilities in spatially developing flows: non-normality and non- 
linearity. Ann. Rev. Fluid Mech. 37, 357-392. 

Cohen, L. S. & Hanratty, T. J. 1965 Generation of waves in the concurrent flow of air and 
a liquid. A.I.Ch.E. Journal 11, 138. 

Craik, A. D. D. 1966 Wind-generated waves in thin liquid films. J. Fluid Mech. 26, 369. 

Delbende, I. & Chomaz, J.-M. 1998 Nonlinear convective/absolute instabilities in parallel 
two-dimensional wakes. Phys. Fluids 10, 2724-2736. 

Delbende, I., Chomaz, J.-M. & Huerre, P. 1998 Absolute/convective instabilities in the 
Batchelor vortex: a numerical study of the linear impulse response. J. Fluid Mech. 355, 
229-254. 

Caster, M. 1962 A note on the relation between temporally-increasing and spatially-increasing 
disturbances in hydrodynamic instability. J. Fluid Mech. 14, 222. 

Gondret, P., Ern, P., Meignin, L. & Rabaud, M. 1999 Experimental evidence of a non- 
linear transition from convective to absolute instability. Phys. Rev. Lett. 82, 1442-1446. 

Hall- Taylor, N. S. & Hewitt, G. F. 1970 Annular Two-Phase Flows. Oxford: Pergamon 
Press. 

Healey, J.J. 2007 Enhancing the absolute instability of a boundary layer by adding a far-away 

plate. J. Fluid Mech. 579, 29. 
Healey, J.J. 2009 Destabilizing effects of confinement on homogeneous mixing layers. J. Fluid 

Mech. 623, 241. 

HOOGENDOORN, C.J. 1959 Gas-liquid flow in horizontal pipes. Chem. Eng. Sci. 9, 205-217. 



36 



L. O Ndraigh, P. D. M. Spelt and S. J. Shaw 



Huerre, P. & Monkewitz, P. A. 1990 Local and global instability in spatially developing 

flows. Ann. Rev. Fluid Mech. 22, 473-537. 
Lecoeur, N., Hale, C. P., Spelt, P. D. M. & Hewitt, G. F. 2010 Visualization of droplet 

entrainment in turbulent stratified pipe flow. In 7th International Conference on Multiphase 

Flow ICMF 2010. 

LlNGWOOD, R. J. 1997 On the application of the Briggs' and steepest-descent methods to a 

boundary-layer flow. Studies in Applied Mathematics 98, 213. 
Mandane, J.M., Gregory, G.A. & Aziz, K. 1974 A flow pattern map for gas-liquid flow in 

horizontal pipes. Intl J. Multiph. Flow 1, 537-553. 
Miesen, R. & Boersma, B. J. 1995 Hydrodynamic stability of a sheared liquid film. J. Fluid 

Mech. 301, 175. 

Miles, J. W. 1957 On the generation of surface waves by shear flows. J. Fluid Mech. 3, 185. 
Monin, A. S. & Yaglom, A. M. 1971 Statistical Fluid Mechanics: Mechanics of Turbulence. 

Cambridge, MA: MIT Press. 
Ozgen, S., Degrez, G. & Sarma, G. S. R. 1998 Two-fluid boundary layer stability. Phys. 

Fluids 10, 2746. 

Pope, S. B. 2000 Turbulent Flows. Cambridge, UK: Cambridge University Press. 

Rees, S. J. & Juniper, M. P. 2010 The effect of confinement on the stability of viscous planar 

jets and wakes. J. Fluid Mech. 656, 309. 
Shampine, L. P., Reichelt, M. W. & Kierzenka, J. A. 1999 Solving Index-I DAEs in 

MATLAB and Simulink. SI AM Review 41, 538. 
O Naraigh, L. & Spelt, P. D. M. 2010 Interfacial instability of turbulent two-phase stratified 

flow: Pressure-driven flow and non-Newtonian layers. J. Non-Newt. Fluid Mech. 165, 489 

- 508. 

6 Naraigh, L., Spelt, P. D. M., Matar, O. K. & Zaki, T. A. 2011a Interfacial instability 
of turbulent two-phase stratified flow: Pressure-driven flow and thin liquid films. Int. J. 
Multiph. Flow 37, 812 - 830. 

6 Naraigh, L., Spelt, P. D. M. & Shaw, S. J. 2012 An analytical connection between tempo- 
ral and spatio-temporal growth rates in linear stability analysis. Eprint: arXiv:1203.1797vl. 

6 Naraigh, L., Spelt, P. D. M. & Zaki, T. A. 20116 Turbulent flow over a liquid layer 
revisited: multi-equation turbulence modelling. J. Fluid Mech. 683, 357 - 394. 

Valluri, P., 6 Naraigh, L., Ding, H. & Spelt, P. D. M. 2010 Linear and nonlinear spatio- 
temporal instability in laminar two-layer flows. J. Fluid Mech. 656, 458-480. 

Ym, C. S. 1967 Instability due to viscosity stratification. J. Fluid Mech. 27, 337. 



